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We study a spin-boson chain that exhibits a local Z 2 symmetry. We investigate the quantum 
phase diagram of the model by means of perturbation theory, mean-field theory and the Density 
Matrix Renormalization Group method. Our calculations show the existence of a first-order phase 
transition in the region where the boson quantum dynamics is slow compared to the spin-spin 
interactions. Our model can be implemented with trapped ion quantum simulators, leading to a 
realization of minimal models showing local gauge invariance and first-order phase transitions. 


I. INTRODUCTION 

Analogical quantum simulators with many-body opti¬ 
cal setups offer us the possibility to replicate the physics 
of condensed matter systems, and also to engineer novel 
exotic quantum phases [1]. In particular, trapped ions [2- 
4] and superconducting circuits [5] are ideally suited to 
implement lattice models of spins coupled to bosons with 
a wide control of spin-boson and spin-spin interactions. 
The resulting family of models that can be directly simu¬ 
lated in those setups include cooperative Jahn-Teller and 
Rabi Lattice Hamiltonians [6-9] . The physical implemen¬ 
tation of those models lead us to the exciting possibility 
to study complex quantum phases governed by the inter¬ 
play between magnetic and vibronic or photonic degrees 
of freedom. Furthermore, the fabrication of arrays of ion 
microtraps open up a new perspective to control lattice 
geometry and particle interactions [10-13]. 

In this work, we introduce a Rabi Lattice model that 
shows a local (gauge) discrete invariance, something that 
takes this model out of the universality classes that we 
typically find in strongly correlated spin-boson lattice 
systems. The implementation of lattice gauge theories 
with trapped ions and superconducting circuits has been 
proposed in recent works [14, 15]. Here we take a dif¬ 
ferent approach to find out the simplest minimal Rabi 
Lattice model that shows local gauge invariance and can 
be implemented in many-body quantum optical setups. 
In fact, spin-boson couplings can lead in a natural way 
to the appearance of a discrete local gauge invariance. 
Consider for example the case of an Ising model, with a 
Hamiltonian of the form 

= -JY. ^Pj+1 ■ ( 1 ) 

J 3 

Local discrete gauge invariance may appear when we 
replace the local field hj by a quantum variable, for ex¬ 
ample, the position operator of a local bosonic field, 

hj g{aj + aj). 


After this substitution, we get an Ising spin model 
where the transverse field is a variable with quantum 
dynamics of its own. This Hamiltonian possess a dis¬ 
crete local symmetry, since it is invariant under a set of 
local transformations defined at each site j, crj ^ 
aj —aj. The model turns out to be a Rabi Lattice, 
where different sites are coupled by an Ising interaction 
between spins. 

This work is organized as follows. Motivated by the 
discussion above, in Section H, we introduce the Ising- 
Rabi Lattice model and its symmetry properties. In Sec¬ 
tion HI we discuss the ground state of the model in some 
limiting cases by using perturbation theory. In Section 
IV we present two variational ansatze to approximately 
find the ground-state of our model: a Born-Oppenheimer 
approximation, valid in the limit in which bosonic de¬ 
grees of freedom are slow compared to the spin dynam¬ 
ics, and a Silbey-Harris approach valid in the limit of 
fast bosonic modes. Those approximations predict a first- 
order phase transition between a pure ferromagnetic Ising 
phase and a dressed ferromagnetic phase of displaced 
bosons. In section V we present numerical results ob¬ 
tained with the Density Matrix Renormalization Group 
(DMRG) method that confirm the validity of the Born- 
Oppenheimer approximation and the existence of a first- 
order phase transition. Section VI presents a proposal to 
implement our model with trapped ions in arrays of mi¬ 
crotraps. Finally we present our conclusions in Section 
VH. 


II. ISING-RABI LATTICE HAMILTONIAN 

We introduce the one-dimensional Ising-Rabi Lattice 
Hamiltonian. Our system consists of N spins arranged in 
a ID chain interacting via a nearest neighbours exchange 
Ising coupling term of strength J (we will assume J > 
0 for definiteness on the following, but the results are 
equivalent if J < 0). Spins are coupled to local bosonic 
modes of energy > 0 by an on-site spin-dependent force 
of magnitude g, 

N N N-1 

-ffiR = (5 X +6' Y G +“i) ~'^Y (2) 

i=i i=i i=i 
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This model possesses a gauge^ (i.e., that acts locally, 
independently on every site) Z 2 symmetry, since it is in¬ 
variant with respect to the transformation prescribed by 


pU) 

' gauge 



SO 


H 

1 ’ gauge 


vj, 


( 3 ) 


that transforms operators ^ 

leaves invariant the Ising coupling term, since cr| ^ crj • 
We expect that this discrete and local symmetry can¬ 
not be spontaneously broken in the ground state of the 
Hamiltonian, a result that is expected to generally apply 
to any local symmetry, and is referred to as Elitzur’s The¬ 
orem [16]. Accordingly, expectation values (aj)GS = 0 
and (crJ)GS = 0 in the whole phase diagram of the model. 

The Hamiltonian (2) also possesses a global Z 2 symme¬ 
try related to the transformation ^ ^2^ whose 

representation in the current space of states is given by 
the unitary operator 


cr? 

V = with V = y] (4) 

i=i ^ 

Since = 0 , the ground state (GS) should ful¬ 

fil (cr|)GS = O 5 unless degeneracy occurs. This global 
symmetry is actually also present in the quantum Ising 
model, where it is spontaneously broken in the ferromag¬ 
netic phase, such that (cr| )gs 7^ 0 in the thermodynami¬ 
cal limit. However, in a finite size quantum Ising chain a 
linear superposition of ferromagnetic states can form the 
ground state, leading to (cr|)GS = 0 foi* finite N. Below 
we show that a remarkable feature of the Ising-Rabi Lat¬ 
tice Hamiltonian is the existence of symmetry breaking 
of the global parity symmetry for finite values of N. 


III. ASYMPTOTIC LIMITS OF THE 
ISING-RABI LATTICE 


to make a explicit choice of basis in the two-fold degener¬ 
ate manifold. To study the stability of the ferromagnetic 
phase, we introduce the spin boson coupling as a pertur¬ 
bation, 

N 

i=i 

and consider its effect upon the degenerate manifold of 
ground states. By applying degenerate perturbation the¬ 
ory, we find that does not lift the degeneracy even at 
finite N (see Appendix (A) for details). This situation 
is in clear contrast with the quantum Ising model, where 
the degeneracy is lifted in the ferromagnetic phase by an 
energy gap scaling like oc [18], with h the value of 
transverse field in Eq. (1). 

By using perturbation theory we calculate the energy 
of any of the degenerate ferromagnetic ground states, in¬ 
cluding the leading corrections induced by the spin-boson 
coupling. 






N-2 2 

(5 + 4J 5 + 2J 


(8) 


Perturbation theory also predicts that states ( 6 ) are a 
good approximation to the ground state of iLjR as long 
as ^ ^ (5 + 4 J. 


B. Dressed ferromagnetic phase 

We consider now the limit g^5 ^ J, where the Ising 
interaction is small compared to the spin-boson cou¬ 
pling and the boson energies. Here we can perform a 
boson-displacement unitary transformation [19] in (2), 
considering as well the rotation x ^ z\ iLjR ^ ^ir = 
with = l/ 2 ^/ 2 (g)f=i(<Tj + <T|), 

and 


A. Ferromagnetic phase 

We discuss the limit 5^ J ^ that leads to a ferro- 
magnetie (F) phase. We define by considering the 
limit = 0 of the IR model. 


N 

^/ = (g)e^^ Sj = ^a]{a]-aj), (9) 

j = l 


SO that the IR Hamiltonian reads iLjR = 
with 5 rp = S ~ Efg‘^/S^ and 


N 

i=i 


N-l 


i=i 


( 5 ) 


^DF = -J 


N-l 

E 

i=i 






V+i 




( 10 ) 


The ground states of consist of the boson vacuum 
and one of the possible ferromagnetic orders (cf. Ising 
model [17]). We will refer to these states as 

N 

|()>F,t) = |0)b(3)l 
i=i 

N 

l<('F,4.) = |0)b(3)l 
i=i 


We have defined operators with Sj ac¬ 

cording to equation (9). The ground states of con¬ 
sist of the vacuum of the bosons in the displaced basis, 
and for any spin configuration. However, this degener¬ 
acy is removed considering the action of the perturbation 
upon these states, 

N-l 

b(0|-ffDF|0)b = y] = f ’ 

i=i 


(6) 
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which shows that the ground states of Hi^i when t ^ 0 
are just the two ferromagnetic states in the x direction. 
Therefore, transforming these states back to the original 
basis, we find the two degenerate ground states of ( 2 ), 

1 ^ 

|0DF,±) = (^(1 - ± (12) 

i=i 

and we will refer to them as dressed-ferromagnetic (DF) 
states. The energy of these states, together with the 
leading order correction induced by the dressing boson 
operators in Hamiltonian (10) is given by 

Euf- - 1 -- J(A^-l)e-4“ -{N-1)—P{a), (13) 

where we have defined 


p{(^) = E 

p=i 


1 e-®""(8Q;2)P 

P p! 


(14) 


Note that P{a) —> 0 if a —> 0, and P{a) —> ( 80 ;^)“^ 
if a oc. The DF state is perturbed by a correction 
{J/6) P{g/6). The latter is negligible if ^ » J in the 
limit g S, and if g^ ^ J5 in the limit g ^ 5. 


C. Qualitative discussion of the quantum phase 
diagram 

The previous considerations allow us to make a conjec¬ 
ture about the phase diagram. We distinguish two cases: 

(i) S J. In this limit, condition g J ensures 

that the F states ( 6 ) are possible ground states of Hiji. 
Following the discussion below Eq. (13), the DF states 
are possible ground states if g ^ VJS. In the interval 
^/J6 < g < the domain of F and DF solutions overlap, 
and we expect a crossover between those energy levels. 
Comparing the F and DF energy, we find that crossover 
at ^ := where we expect the appearance of a 

first order F-DF transition. 

(ii) 6 ^ J. Here, F states are valid ground states if 
g 6, whereas DF states are valid ground states for 
any value of g, as follows from the discussion below Eq. 
(13). In the interval ^ ^ E and DE solutions overlap, 
however, here the DE state continuously converges to the 
E state. Thus we expect a continuous transition from the 
DE to the E solution. 

Putting together all previous arguments, we expect 
that presents a first order quantum phase transition 
along the critical line gc{S,J)^ featuring a jump from the 
E to the DE ground states in the regime of low boson 
energies J —> 0. This is in clear contrast with the quan¬ 
tum Ising chain with a transverse field, where there is no 
coexistence of the ferro- and paramagnetic phases at nei¬ 
ther side of the (second order) phase transition. In the 
however, there is a coexistence of the phases already 
addressed if \/^ ^ g J. Eurthermore, this last set of 


inequalities cannot be longer fulfilled if5':^J^ and there¬ 
fore the discontinuous behaviour is bound to disappear 
for a given S ^ J. We have summarized these considera¬ 
tions in Eig. 1, where we choose as order parameter the 
average boson number 


«=;^E(aJai) (15) 

i=i 

to capture the sudden change from the boson vacuum 
state (E phase) to a displaced state (DE phase). 



FIG. 1. Scheme of the phase map depicting the disappearance 
(at the dot) of the discontinuous jump in the number of bosons 
along the critical line (solid) for a given value of S, g. The 
dashed line represents no boundary but a continuous tran¬ 
sition from the ferromagnetic to the dressed-ferromagnetic 
phase. 


Below we consider two different mean-field descriptions 
that give some physical insight on the phases and the 
transition of the problem. In addition, they will be vali¬ 
dated afterwards by a exact numerical DMRG diagonal- 
ization. 


IV. VARIATIONAL METHODS 


A. Born-Oppenheimer approximation ((5 < J) 


The classical limit of the model is attained in the 
regime of very high number of bosonic excitations, which 
is expected at J —> 0. In this limit, ladder operators can 
be treated as classical variables, aj aj G C, and Hi^ 
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is reduced to a spin Hamiltonian, 

Hbo = 

N N N-1 

^ K + “i) crpi+i-( 16 ) 

j = l j = l j = l 

Eq. (16) describes an Ising chain in a transverse field, 
for which an exact ground state |4^i(aj)) can be found 
[17]. Without loss of generality we can assume aj to be 
real. We devise a variational ansatz by calculating the 
mean value of Hbo^ whose ground state energy can be 
written as EBo{{aj}) = where 

^ 1 , 0 is the ground state energy of the quantum 
Ising chain (1) with transverse fields hj = 2gaj and inter¬ 
action strength J. The corresponding variational wave- 
function is 


N 

I^'bo) = l^i(ai)) \otj)- (17) 

i=i 

This method is a self-consistent approach that resem¬ 
bles the Born-Oppenheimer approximation in Molecular 
Physics [20]. In that context, the degrees of freedom of 
the positions of the nuclei enter the electronic Hamilto¬ 
nian as parameters in the same way the boson ampli¬ 
tudes appear in the spin Hamiltonian (16). We notice 
that, due to the underlying gauge symmetry in the Hib 
Hamiltonian, a variational solution of the form (17) can 
be transformed into a solution with the same energy if 
we change locally the sign of the displacement and 
simultaneously transform crj ^ • There are thus 2^ 

possible solutions, given by the values aj = Sj\aj\^ with 
Sj = ±1. 

In order to make best use of the analytical results for 
the solution of (16), we assume N ^ oc and aj —a 
(in the thermodynamic limit the system is homogeneous, 
whereas the minus sign is chosen for analytical conve¬ 
nience), so the energy Ebo of the ground state of Hbo 
is 


Ebo 


Sa^-2ag-{lEX)E 

TT 


4A 


(i+A)2 




where E is the complete elliptic integral of the second 
kind. We are interested in the value of the parameter a 
for which the energy attains a minimum; we will refer to 
this point as o^o, and its value together with the exact 
solution of the spin problem will define the mean-field 
ground state. 

A quick inspection of (18) reveals that as a function of 
a, the energy is minimum exclusively at the origin unless 
there are values of J^g that shift its position to a finite 
value a 0 (see Fig. 2). Bearing this in mind we carry 
out the Taylor expansion of the energy around a = 0, 
that leads to 


^=-j+(S-^)a^ + 0{a^), (19) 



a 


FIG. 2. Mean field energy (18) for (5 = J = 1 as a function 
of the parameter a and different values of the spin-phonon 
coupling g. Note that close to the origin there is a curvature 
change for a given g > gc- This point marks the criticality 
condition. 


and predicts a minimum for o^o 7 ^ 0 (o^o = 0 ) whenever 
g^ > SJ := g^ {g < gc). We interpret that the system 
undergoes a phase transition at that point: below the 
critical line bosonic excitations are inhibited (q^o = 0 ), 
and the spins point in the Ez or — 2 : directions; above 
the ground state changes abruptly to allow an arbitrary 
number of bosons, whereas the spins point in the direc¬ 
tion determined by the Ising ground state for a trans¬ 
verse field of magnitude 2^05^. Furthermore, in this latter 
regime we can estimate the value of assuming g ^ J 
-for fixed a, S-^ which gives 




( 20 ) 


From the previous discussion we can extract the or¬ 
der parameter n = q^q, that we shall compare with the 
DMRG results to assess the validity of the previous ap¬ 
proximations. 

The Born-Oppenheimer solution converges to the DF 
states in the limit g ^ S. In this limit one can easily 
show that the optimal values are \aj \ = g/S. We can re¬ 
store the Z 2 gauge symmetry by considering a symmetric 
superposition. 


N 

i‘tB’'o) = 5W5 E l*■(»4)>(8)l•4)■ 

7 = 1 

S,=±l 

such that we recover the solution |0 df,+)- The solu¬ 
tion I 0 DF,-) would correspond to the antisymmetric lin¬ 
ear combination of the former states. 


N 
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B. Silbey-Harris-type ansatz (6 > J) 

In order to investigate the continuous transition regime 
mentioned in Fig. 1, we are going to consider a displaced 
trial wave-function whose distance away from the origin 
in phase space is no longer fixed, rather the variational 
parameter [21]. This approach has been recently shown 
to yield an accurate description of the quantum phase 
diagram in Rabi Lattice models [9, 22]. 

Specifically, we take the IR Hamiltonian in the rotated 
basis X ^ and compute its energy upon the wave- 
function 

N N 

I^sh) = |0)b (g) It.),-, S{7j) = 7?| ^ a;(at - a), 

i=i i=i 

( 22 ) 

where the parameter 77 continuously interpolates the dis¬ 
placed solution between 0 and g/6 for fixed values of 
these. The Silbey-Harris energy reads 

^sh(^) = (r/" - 2r?) - J{N - , (23) 

which along with the condition dEsu/dr] = 0 for a given 
77 = 770 , leads to the optimal value for the order parameter 
n = {r]og/S)‘^ within this framework. 

This ansatz resembles the exact IR Hamiltonian solu¬ 
tion if (5 ^ J, because in that case the ground state is 
one of the dressed-ferromagnetic eigenvectors (cf. section 
HIB). However, it turns out that it also predicts a first 
order phase transition when extrapolated to the 6 J 
regime. This supports the fact that exhibits a sud¬ 
den ground state change in this latter case (see Fig. 3). 

We present the predictions of the Silbey-Harris solu¬ 
tion, focusing on the fact that the discontinuity of n at 
the transition disappears between the regimes 6 J and 
6<^J. 


V. DMRG RESULTS 

In this section we present quasi-exact numerical calcu¬ 
lation of the ground state properties of the IR Hamilto¬ 
nian for a chain of N = 50 spins, obtained by means of 
the DMRG algorithm [23]. Some remarks are in order 
before proceeding to the results. First, we have to intro¬ 
duce a cut-off, Nc^ in the maximum Fock state of local 
bosonic modes in the DMRG algorithm. This imposes 
some limitations in the description of the DF phase in 
the limit 6 g. Here, due to the low energy cost of 
bosonic excitations, the ground state wavefunction has 
non-negligible projections upon many different occupa¬ 
tion states. Thus, an accurate description may require 
high values of Nc that are beyond our computational ca¬ 
pabilities. In this work we use Nc = 10, and present 
exclusively DMRG results fulfilling 2n < Nc. Second, 
we stress that finite-size effects in our calculations lead 
to a smearing of discontinuities at the first order phase 



FIG. 3. Silbey-Harris mean boson number n for different 
values of (5, J = 1 and iV = 50 sites. 


transition, which strictly speaking takes place in the ther¬ 
modynamic limit only. 

However, our finite-size results are consistent with the 
occurrence of a first order phase transition in the regime 
of slow boson dynamics. To assess this phenomenology, 
we focus on the behaviour of the mean boson number n. 
As depicted in Fig. 4, n shows a sudden change when 6 
lies deep in the regime 6 whereas the discontinuity 
vanishes for S > J (we set units such that J = 1). To 
quantify better that discontinuity and to place accurately 
the position of the phase transition, we have computed 
the numerical derivative of ti as a function of g, see Fig. 
5. We observe that for S below J the numerical derivative 
inversely scales with 5. This results is consistent with the 
sudden change in the ground state between the n O^F 
phase, to the displaced vacuum of the DF phase, where 
n = aQ S~‘^ according to Eq. (20). Increasing the 
values of 6 leads to a disappearance of any peak in the 
numerical derivative. We conclude then that the discon¬ 
tinuous behaviour is only unveiled in the limit 6 J, 
because any signature is lost when 6 > J. 

We have also compared the exact results with the vari¬ 
ational approaches. Let us start by checking the accuracy 
of the Born-Oppenheimer approximation, which works in 
the limit 6 J^g. To this end, we look for a closer re¬ 
semblance between the BO solution and the exact diag- 
onalization for decreasing values of 6 (cf. Fig. 4). Ac¬ 
cordingly, we see that the smaller the S, the nearer the 
BO prediction for the number of bosons n lies to the 
DMRG observable. This is also true in the case of the 
derivative of the number of bosons where, in contrast to 
the Silbey-Harris ansatz, the BO approximation quanti¬ 
tatively predicts the height of the derivative when J —> 0. 
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FIG. 4. Mean boson number prediction for Born- 
Oppenheimer (dashed lines), Silbey-Harris-type ansatz 
(dashed-dotted lines) and DMRG (solid lines) of a A/" = 50 
sites chain, and J — 1. For the DMRG method we set a renor¬ 
malization dimension D = 10, on-site boson cut-off A^c = 10 
and local dimension d — 2- (we follow the notation of [24] ). 


Regarding the Silbey-Harris approach, Fig. 4 shows 
that it correctly describes the existence of the disconti¬ 
nuity. However, this solution must also give a suitable 
description of the phase with (5 ^ J, as we know that 
the dressed-ferromagnetic phase consists of a displaced 
state. We have therefore run simulations for bigger val¬ 
ues of 5 and g (cf. Fig. 6) and compared them with 
the SH ansatz, that effectively coincides with the exact 
solution when 5^g ^ J. 

In Fig. 7 we present the scaling of the critical line 
with the parameter S, in the regime 6 < J. It has been 
obtained from the position of the maxima of the deriva¬ 
tives of n as a function of g, for different values of S. This 
allows us to calculate the function gdS)^ defining the crit¬ 
ical line. Our results yield a power law, e.g., gc ^ for 
fixed J, with the exponent a = 0.66. The quasi-exact nu¬ 
merical result departs from the BO approximation, that 
predicts gc = \/^, that is, a = 1/2. 

We have also studied signatures of the first order phase 
transition in the correlation length. Let us define the spin 
correlation functions, 

= (24) 

Our calculations show that Cz{i^j) oc along 

the whole phase diagram, where x is the correlation 
length. The exponential decay is observed even close 
to the first order phase transition in the regime 6 < J. 



FIG. 5. Derivatives of the average boson number for different 
values of S {N = 50, J = 1). Note that the DMRG diago- 
nalization (solid lines) gets closer to the Born-Oppenheimer 
prediction (dashed lines) for decreasing S, whereas the Silbey- 
Harris ansatz (dashed-dotted lines) improves for bigger values 
of the boson energy. The step for the derivatives in all cases 
is the same and stems from the precision used in the DMRG 
diagonalization: Ag = 0.02 • J. 


This is consistent with our picture of the transition as a 
level crossing: F and DF states are both close to eigen¬ 
states of iLiR at the critical point, and both of them show 
exponentially decaying correlations. This is in clear con¬ 
trast with what one would expect in a second order phase 
transition [25]. On the critical line, S can be identified as 
the energy gap separating the ground state sector from 
the lowest energy excitations. We thus expect that the 
correlation length on the critical line, Xc, must be a de¬ 
creasing function of S. Our DMRG calculations confirm 
this picture (Fig. 8), and yield the scaling Xc oc 1/5 (Fig. 
9). 


VI. IMPLEMENTATION OF THE ISING-RABI 
LATTICE HAMILTONIAN WITH TRAPPED 
IONS 

In this section we discuss an eventual realization of the 
IR Hamiltonian in state-of-the-art trapped ion set-ups, 
where highly accurate state preparation and readout is 
currently achievable [3]. In these systems, two electronic 
levels of every ion are chosen to be regarded as the spin 
degrees of freedom, whereas the quantized oscillations of 
the ions (phonons) give rise to the bosons. Then, spins 
and phonons are coupled through optical forces. We de¬ 
vise using a linear array of microtraps [10-13] instead of 
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FIG. 6. DMRG average boson number (solid lines with 
symbols) vs. Silbey-Harris ansatz (dashed line) for > J 
(iV = 50, 1). 



log 5 


FIG. 7. Linear fit of the critical line gc(S) from the DMRG 
results {N = 50, J = 1). We depict both the logarithm of 
S,g. The result is coherent with a power law decay gc ^ S^, 
with a ~ 0.66. 


the more usual Paul traps [26] . Individual traps are spe¬ 
cially suitable for our purpose, as their frequencies can 



FIG. 8. Gorrelation lengths y obtained as the slope of 
— Cz(25, 25 + j)/N), across the critical line as a 

function of g, for different values of (5 (A^ = 50, J = 1). 



FIG. 9. Fitting of Xc^ fo a line, from the DMRG with N = 50, 
J — 1. The results are coherent with the fact that the gap 
A scales linearly with the bosons energy S along the critical 
line. 


be independently tuned, and the motion of every ion can 
be made resonant with a different laser force. 

There exist as well other experimental set-ups for the 
















simulation of the Hamiltonian such as supercon¬ 

ducting qubits [5] or Rydberg atoms [27]. For example, 
in this latter case the spin-spin interaction -between the 
two level systems made up of the ground and (very high) 
excited state of every atom in the sample- is directly 
induced by the electromagnetic interaction between the 
electronic states. The spin-boson coupling can be in¬ 
troduced by the action of lasers as in the trapped ions 
experiments. 

A. Description of the set-up 

We discuss the set-up in three parts, consisting of the 
three different terms in i^iR. The reader is referred to 
[3, 26] for further details of the implementation. 

1. Phonon Hamiltonian 


We need to get rid of the hopping terms cil,jCia,i in 
^phonon, at Icast for one direction a, in order to give 
rise to the local boson contribution in the IR Hamil¬ 
tonian, e.g., Let us choose for this purpose 

the transversal modes along the x axis. Then, we pro¬ 
pose using different trap frequencies uJxj to make hop¬ 
ping events in (27) fast rotating compared to the on-site 
energies. Specifically, if ions j and I are subjected to fre¬ 
quencies ujxj and (jJx,i, the terms would rotate 

with e-xjp[—it{(jJx,j — ^x.i)] in the interaction picture for 
the motion. The so-called Rotating Wave Approximation 
(RWA) prescribes that such terms are negligible as long 
8iS tj I ^ \^^x,j —^x,i\‘ Assuming this is the case, hopping 
terms in iLphonon can be safely ignored, and we are led 
to the term 

N 

Hx — 

i=i 


Let us assume a linear array of traps forming an ion 
chain along the z axis. Furthermore, we consider a con¬ 
stant separation do between traps. The (quantum) posi¬ 
tion of the ions can be written as 


Tj — 6txjx Svyjy -|- (^Zj + 5tzj^z^ (2^) 

where operators Svaj stand for the displacements off 
their equilibrium coordinates (0, 0,2;^). Motion along the 
y axis is not relevant for the simulation, and will be omit¬ 
ted in the following. Displacements between different 
directions are decoupled assuming effectively harmonic 
trapping potentials, and approximating the Coulomb in¬ 
teraction up to second order in Therefore, ions are 

subjected to the effective potential 






C(yC 


72 


a 


■yC 13 


{SVaJ - Sra,lf 


(26) 


In this expression Cx = I, Cz = —2, m is the ion mass, 
e the electron charge in CGS units, and are trap 
dependent frequencies. The corresponding Hamiltonian 
can be canonically quantized expressing the positions and 
momenta in terms of creation and annihilation operators, 
so that [19, 28, 29], (we take H = 1) 


we were aiming for. The common frequency ujxj ^ 
can be achieved by means of local laser detunings, dis¬ 
cussed later on. The motional coupling between different 
traps decays fast as a function of the ion-ion distance, 
tji ^ l/\Zj — zf\^. Thus, it is only necessary to eliminate 
the coupling between nearest or next-to-nearest neigh¬ 
bour ions, since longer-range terms will give negligible 
contributions. 

Regarding the motion in the 2 ) direction, we set oozj 
ujz> Since trap frequencies along x and 2 ) are indepen¬ 
dently and locally tunable, this choice can be made at 
no expense of the previous discussion. The Hamiltonian 
(27) reads then 


AT-l 

Hz — ^ ^ (^0) 

n=0 


in the basis of collective modes of motion az^n = 
with normal frequencies uJz,n- This term 
does not occur in TLir as we aim at a regime where 
{al ^az^n) — 0- Nonetheless, their (virtual) exchange cre¬ 
ates the spin-spin interaction [19]. 


2. Spin-spin interaction 




phonon 


where 


« 3 3 1^3 

(27) 


= Z 




2m(w„jWa,;)l/2|^0 _ 2;0|3 ' 


(28) 


In (27) is already assumed that ^ such 

that corrections to on-site frequencies stemming from the 
dipolar interaction, or phonon non-conserving terms, are 
negligible. 


Implementing the exchange term of TLir relies on in¬ 
ducing a spin-spin effective coupling. Let us assume a 
laser field, lying along the direction of the linear ar¬ 
ray of traps, with momentum /\kz and frequency = 
— ^z,n- Here 6z,n stands for the laser detuning from 
the n axial normal mode. Differential a.c. Stark shifts 
stemming from the off-resonant components of the atom- 
light interaction give rise to a spin-dependent a^-force 
[30] of the form 

-f^z—force ~ 9z ^ ^ <^7 (y^j,n^z,n T H.C.^ , (^ 1 ) 

3,n 
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FIG. 10. Scheme representing trapped ions in linear array of 
microtraps as electrodes printed over a surface. Solid arrows 
represent the laser fields acting on the ion chain. We indicate 
the trap frequencies at every site. 


with coupling strength Qz. This Hamiltonian is time in¬ 
dependent because we have moved to a rotating frame, 
where phonon frequencies are shifted, ujz,n = 

5z^n- Performing a transformation to a displaced basis 
(see e.g. [19]), the previous force takes the form of the 
effective spin-spin interaction 

-f^exchange = ^ ~ "P TP’ 

3,1 ~ ' 

for suitable detunings and shapes of the axial modes spec¬ 
trum. This interaction acts effectively as a first neigh¬ 
bours ferromagnetic coupling of magnitude J. 


3. Local spin-phonon coupling 

Local spin-phonon couplings in iLjp require driving si¬ 
multaneously red and blue sideband transitions [26] for 
the transversal oscillations. However, as we have already 
discussed, ujx,j are different among close traps. This 
means that matching the resonance conditions for the 
spin-dependent forces requires as many laser wavelengths 
as different trapping frequencies. Let us consider the ar¬ 
ray of traps as consisting of N/n^ n G N sequential sets 
of traps. Within these, neighbouring traps frequencies 
are different. We set a constant difference between one 
trap and the next, ujx,j — All the sets 

have the same arrangement of n frequencies, and they 
appear one after the other along the chain. Let us call 
these frequencies uJx.h • • • ^^x,n- Any frequency can be 
written then as 0Jx\j]^ where [j] = {j — 1) mod n + 1. 
Now, we apply n laser fields transversally to the chain, 
with mutual detunings Aux ,..., (n — l)Aujx> Because 
of this frequency difference, they can address the whole 
chain at the same time. In this way, the matching condi¬ 
tion only happens between a given laser with, let us say 


[j] = ^0 + ^x,[j] ~ ^x,[j]^ arid the ions that are trapped 
at frequencies ujx,[j] (<^o is the spin transition frequency). 
This gives rise to the cr^-force 

N 

-H'x-force(i) =9'^ O-J (33) 

i=i 

where g = iftx,[j]Vx,[j]^ fhe laser Rabi frequency and 
Lamb-Dicke parameters of the coupling, respectively. We 
rely on the local dependence of ^x,[j] fo achieve a homo¬ 
geneous g along the chain, as r]x,[j] depend on the on-site 
trap frequencies. 

Finally, moving into a rotating frame with frequencies 

i?x-force(i) 

-force (0), SO that 

N 

-f^x—force ~ 9 ^ ^ {^x,j (^d) 

i=i 

Since laser detunings are site-dependent, they can be 
shifted to give common on-site phonon energies (5, Vj, 
which leads to 


N 

Hx = (35) 

i=i 

as the effective phonon energy contribution. 

The IR Hamiltonian is eventually implemented as the 

sum of Hx^ -^exchange S-Ild Hx—force- 

B. Trapped ions experimental parameters 

We consider the traps separated by a distance do = 
30 /im, every of them containing one ^Be+ ion. We es¬ 
timate \z^- — Zi \ = do in (28). Eq. (31) holds only if 
m.?cyiiriz,n = Akz/y/2muJz,n=o ^ 1- We propose a com¬ 
mon ujz = 500 (27r) KHz for all traps, which leads to 
— 29 (27r) KHz and to uJz,n=o — 431 (27r) KHz 
for the ground state COM frequency of the axial modes 
band. Therefore, a laser wavelength 870 nm would 
give r]z,n=o — 0-26 for beams on axis with the traps. 

The magnitude of the exchange in (32) is J 
tjj-\-i9z/^z,n=o^ where gz has typical values 100 (27r) KHz 
[3], whereas we impose 6z,n=o — in order to neglect 
residual spin-phonon couplings [19]. This renders the 
value J 'Zi 1 (27r) KHz, which is the lowest energy scale 
involved in the simulation. 

The number of different frequencies ujx^j fixes an er¬ 
ror bound for the simulation. Ions trapped at equal 
frequencies are coupled by a residual dipolar interac¬ 
tion, whose magnitude is maxj(tj^^-^^). We aim at mak¬ 
ing it very small with respect to the rest of parameters 
in idiR. Then, processes with energies ^ maxj(t^j^^) 
will be systematically neglected. For the sake of con¬ 
creteness, we address the example of n = 3. Assum¬ 
ing uJx\i] = 10 (27r) MHz, uJx\ 2 ] = 9 (27r) MHz, and 
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^x,[3] = 8 (27r) MHz, we have 0.9 (27r) 

KHz. This amount scales with the distance, so that 
^ ^ 33 (27r) Hz. Accord¬ 
ingly, we prescribe d^g/j ^ as the condi¬ 

tion to be fulfilled to safely neglect residual couplings. 
Furthermore, with the former choice of parameters, the 
RWA condition is also fulfilled, as — 

>^x,i\) - 10“^, ^ = j + 1, • • • ,n. 

Regarding the spin-boson interaction, we consider laser 
beams with effective wavelength Ag — 320 nm acting 
transversely to the traps’ axis. Thus, the Lamb-Dicke 
parameters are max/y^^^j 0.16. Typical values for g are 
again of the order of 100 (27r) KHz. The energy of the 
transverse phonons is set by locally detuning from ujx,j 
to the common value 5 for every site, and it can chosen 
such that 5 ^ g^ as we have theoretically studied. 

In order to probe the phase transition we propose 
preparing the ferromagnetic phase by cooling to the 
ground state of the phonons, while optical pumping to 
the \iz)j spin state, where \iz)j is one of the qubit 
states. An adiabatic protocol crossing the critical line 
would require evolution times of the order of the in¬ 
verse of the smallest of the parameters, which lies around 
^ 23 /is. 


VII. CONCLUSIONS 

We have introduced the Ising Rabi Lattice model, 
that consists of a generalization of the single particle 
Rabi model that includes Ising couplings between spins. 
Our model departs from the Ising universality class, and 
presents a discrete gauge symmetry. We have used sev¬ 
eral approximations and perturbative arguments that 
predict a quantum phase diagram divided in two parts: 
(i) Slow boson regime {6 J) in which a first order 
phase transition separates a ferromagnetic phase from a 
phase with a dressed-ferromagnetic phase, (ii) Fast bo¬ 
son regime {6 J), were the transition between the F 

and DF phases is continuous. This picture is consistent 
with quasi-exact numerical calculations with the DMRG 
method. Our model can be implemented with trapped 
ions in arrays of microtraps, leading to the implementa¬ 
tion of gauge symmetries and first order phase transitions 
in this system. 


ACKNOWLEDGEMENTS 

P.N. acknowledges H. Takahashi for its help in details 
of the implementation. D.P. acknowleges J.J. Garcfa- 
Ripoll for suggesting the application of the Silbey-Harris 
variational approach in this work. Work supported by 
the EU Marie Gurie Gareer Integration Grant 630955 
NewFQS. 

Appendix A: Survival of the two-fold degeneracy up 
to any finite order in perturbation theory 

In the ^ = 0 limit of i^iR, its ground states are those 
of 

N 

i=i 

which fulfil 

N N I 

|0)b(g)|t.),|0)b(g)|;.) b (A2) 

i=i i=i J 

where 1. c. stands for every independent linear combina¬ 
tion of these vectors. Because of the two-fold degener¬ 
acy, we are allowed to choose as ground states any two 
elements of (A2), so for the sake of simplicity we will 
consider that |(/)f) is any of the ferromagnetic orders 

N N 

\h) ■■= |o)b (g) It.) , \h) ■■= |o)b (g) It.) . (A3) 

i=i i=i 

If we consider now the perturbation 

N 

(A4) 

i=i 

the first step in building corrections to |0 f) consists in 
finding the correct linear combination, i.e., the weights 
in |(/)f) {g ^ 0)) = "that continu¬ 

ously matches the ground state for ^ = 0. This is accom¬ 
plished by means of degenerate perturbation theory (see 
[31] for an account of the Brillouin — Wigner approach), 
which studies the effect of the perturbation within the 
degenerate subspace (A3): if is such that it lifts the 
degeneracy at a given order n, this procedure provides 
two new eigenvectors whose energies are the eigenvalues 
of the following secular equations 


N-l 

i=i 


(Al) 


E. 






GS 


c-|- — ^ ^l'?^t)ct + ('/’tiFE 


H' 




Hi, 


RgsHj 
H' 




RgsH, 




(A5) 
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In the former expressions are the ground 

state energies for ^ = 0 and i?GS = ~ ~ 

|0f)(0f|) is known as the resolvent. 

In the present case, we are going to show that the de¬ 
generacy is not lifted at any finite order in perturbation 
theory. We note that equations (A5) give only one solu¬ 
tion for Eq^ if and only if both 


{4>t\Y 


H' 


ihlYZ 


■ -Rcs-ffi 
HL 


tIH) = ih\Y. 


_Hf _u \ 

D TJf 1^1/’ 

-ttGS-^F 


(A6) 


RgsH, 




do hold. The first of these conditions is trivially fulfilled 
because of parity arguments, but the second must be 
computed explicitly. It turns out that it holds as well, 
because all the matrix elements in 


Hi. 


H' {RgsH^T |<^t) 


/ \k I 


/c=0 


(A7) 


are zero. To show this, let us write the generic form of a 
n-th order contribution to the previous sum (denomina¬ 
tors can be neglected for this discussion), 

N 

5” (g) b(0|(a] + |0)b U.l |t.), (AS) 

i=i 

with We note that the boson displacement 

terms are diagonal only in the event of even values of ev¬ 
ery Uj. However, only an odd value of all the nj would 
give a non-zero contribution from the spin part, because 
in another case the tunneling matrices are equal to the 
unit matrix. Since both contributions cannot be simulta¬ 
neously different from zero, we conclude that the second 
condition in (A6) is also fulfilled. 

According to these previous considerations, the two¬ 
fold degeneracy in the ground state is not lifted in any 
finite order of perturbation theory, which in turns trans¬ 
lates into the fact that degeneracy remains for any finite 
value of g [25]. Therefore, perturbative corrections must 
be carried upon any element of (A3) by means of con¬ 
ventional non-degenerate perturbation theory. 
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